Biostar：课程25、26

查看原文

其他

Biostar：课程25、26

原创 2017-09-09 István Albert 生信媛

翻译：生物女博士

注：本系列课程翻译已获得原作者授权。

# 为作者的注释

#* 为译者的注释

Lecture 25 - Bedtools简介

# 安装bedtools

#* 后台有人发消息来问，你怎么知道你要找的这个软件的链接在哪里？

#* 其实，很多官网上，都是有的。跟你下个office、PS等这些软件是一样的道理。

#* 比如今天要学习的bedtools

#* 首先，百度或者google关键词：bedtools

#* 搜索得到的词条，前两个都是官网的（有时候不一定排在很前面）。

#* 进入官网以后，就找到了。可以看到是在github上存放的。

#* 先大致看一下网页，对于github，可以点击“clone or download”点击下载或者获得下载的链接。

#* 具体下载方式见下图：

#* 在本软件的下载中，在页面下方，还有另外一个安装包下载链接：

#* 留给有兴趣的读者自行探索。

#* 其他也是相类似的搜索。如果是R包，则可以到CRAN（https://cran.r-project.org/）上的Packages（https://cran.r-project.org/web/packages/index.html）或者是bioconductor（http://www.bioconductor.org/）网站去找。

#* 好了，进入今天的正题。

# 安装bedtools

cd ~/src
curl -OL https://github.com/arq5x/bedtools2/releases/download/v2.22.0/bedtools-2.22.0.tar.gz
tar zxvf bedtools-2.22.0.tar.gzcd bedtools2
make
ln -sf ~/src/bedtools2/bin/bedtools ~/bin/bedtools

# 演示版的bed文件 (demo.bed)

KM034562    100    200    one    0    +
KM034562    400    500    two    0    -

# 我们的基因组文件（genome.txt）

KM034562 18959

# 没有链特异性的运算

bedtools slop -i demo.bed -g genome.txt -l 10 -r 0

# 有链特异性的运算

bedtools slop -i demo.bed -g genome.txt -l 10 -r 0 -s

# 把BED转成对应的GFF

# 这并非是真的正确地把BED转成GFF。

cat demo.bed | bioawk -c bed '{print $chrom, ".", ".", $start+1, $end, $score, $strand, ".", "." }' > demo.gff

# 看！它与其他格式可以很好地协同工作！那怎么可能呢？好神奇！

bedtools slop -i demo.gff -g genome.txt -l 10 -r 0 -s

# 两侧的运算。把结果重定向到一个文件中。

bedtools flank -i demo.gff -g genome.txt -l 10 -r 0 -s > flank.gff

# 填充运算。

bedtools complement -i demo.gff -g genome.txt > complement.gff

# 获取埃博拉基因组。

#* KM034562.gb文件在之前的课程其实已经有给过大家了。

#* 或者也可在这里下载，链接：https://share.weiyun.com/e880cca38636c8d90780fb487050bf3d （密码：HhCUk9）

efetch -id KM034562 -format gb -db nucleotide > KM034562.gb
readseq -format=FASTA -o ~/refs/852/KM034562.fa KM034562.gb
readseq -format=GFF -o KM034562.gff KM034562.gb

# 看一下提示，这个工具可以用于哪些格式。

bioawk -c help

# 从整个文件中提取基因。

cat KM034562.gff | bioawk -c gff ' $feature=="gene" { print $0 } ' > genes.gff

# 序列提取。获取与间隔相对应的序列。

bedtools flank -i genes.gff -g genome.txt -l 10 -r 0 -s > flank.gff

#* 在这一步，使用最新版本的bedtools会有奇怪的报错。而用作者提供的这个旧版本则不会。

# 提取与genes.gff的间隔相对应的序列。

bedtools getfasta -fi ~/refs/852/KM034562.fa -bed genes.gff -fo genes.fa

# 来看一下结果文件。

head genes.fa

Lecture 26 - Bedtools指南

# 来到ecture 23，我们进行了埃博拉的短reads比对

cd ~/edu/lec23

# 创建一个间隔文件，叫region.bed

# KM034562 1000 2000

# 我们可以用这个间隔文件去分割数据。

bedtools intersect -a bam/SRR1553593.bam -b region.bed > region.bam
samtools index region.bam

# 我们也可以把数据作为bed文件输出。

bedtools intersect -a bam/SRR1553593.bam -b region.bed -bed > region.bed

# 又或者，我们可以把bam文件转成bed文件。

bedtools bamtobed -i bam/SRR1553595.bam > SRR1553595.bed

# 我们将跟着bedtools的指南来做

http://quinlanlab.org/tutorials/cshl2014/bedtools.html
mkdir ~/edu/lec26
cd ~/edu/lec26

# 获取所有的数据，并解包。

curl -O http://quinlanlab.cs.virginia.edu/cshl2013/maurano.dnaseI.tgz
curl -O http://quinlanlab.cs.virginia.edu/cshl2013/cpg.bed
curl -O http://quinlanlab.cs.virginia.edu/cshl2013/exons.bed
curl -O http://quinlanlab.cs.virginia.edu/cshl2013/gwas.bed
curl -O http://quinlanlab.cs.virginia.edu/cshl2013/genome.txt
curl -O http://quinlanlab.cs.virginia.edu/cshl2013/hesc.chromHmm.bed

# 这个数据集含有更多的数据。

# 胎儿的组织样品，包括脑、心脏、肠道、肾脏、肺、肌肉、皮肤以及胃

tar -zxvf maurano.dnaseI.tgz

# 我们将在课堂上跟着指南做一点实际操作

# 并给作业留一些建议

# 从注释文件中，选取启动子。

cat hesc.chromHmm.bed | grep Promoter > promoters.bed

# 找到跟每个exon最近的启动子。

# -d 表示打印列

#* 多的一列数值是-a 和 -b 两者最近的距离。（见下图）

bedtools closest -a exons.bed -b promoters.bed -d | head -1

# 找到覆盖了最多外显子的CPG岛